A framework for estimating and visualising excess mortality during the COVID-19 pandemic

COVID-19 related deaths underestimate the pandemic burden on mortality because they suffer from completeness and accuracy issues. Excess mortality is a popular alternative, as it compares observed with expected deaths based on the assumption that the pandemic did not occur. Expected deaths had the pandemic not occurred depend on population trends, temperature, and spatio-temporal patterns. In addition to this, high geographical resolution is required to examine within country trends and the effectiveness of the different public health policies. In this tutorial, we propose a framework using R to estimate and visualise excess mortality at high geographical resolution. We show a case study estimating excess deaths during 2020 in Italy. The proposed framework is fast to implement and allows combining different models and presenting the results in any age, sex, spatial and temporal aggregation desired. This makes it particularly powerful and appealing for online monitoring of the pandemic burden and timely policy making.


Introduction
Estimating the burden of the COVID-19 pandemic on mortality is an important challenge (Weinberger et al., 2020). COVID-19 related deaths are subject to testing capacity and changes in definition and reporting, raising accuracy and completeness considerations (Aburto et al., 2021;Konstantinoudis et al., 2022). In addition, COVID-19 related deaths give an incomplete picture of the burden of the COVID-19 pandemic on mortality, as they do no account for the indirect pandemic effects due to, for instance, disruption to health services (Kaczorowski and Del Grande, 2021). Alternatively, excess mortality has been extensively used to evaluate the impact of the COVID-19 pandemic on mortality (Rossen et al., 2020;Islam et al., 2021;Kontis et al., 2020;Konstantinoudis et al., 2022;Verbeeck et al., 2021).
Excess mortality is estimated by comparing the observed number of deaths in a particular time period with the expected number of deaths under the counterfactual scenario that the pandemic did not occur. Typically this counterfactual scenario is calculated using a comparison period, for instance several previous years (https://www.euromomo.eu/). Calculating accurately the expected deaths requires accounting for factors such as population trends, seasonality, temperature, public holidays and spatio-temporal dependencies. While accounting for the above-mentioned factors, most studies to date have estimated excess mortality at national level (e.g., Rossen et al., 2020;Weinberger et al., 2020) and a few have looked across different countries reporting important geographical discrepancies (Islam et al., 2021;Kontis et al., 2020Kontis et al., , 2021. While national level estimates of excess mortality give valuable insights into the total number of excess deaths, they do not allow the evaluation of within country geographical disparities. Such information is essential to understand the country's transmission patterns and effectiveness of local policies and measures to contain the pandemic (Kontopantelis et al., 2021). Temporal trends in excess can substantially differ across regions of the same country, which makes national-based comparisons even more challenging (Blangiardo et al., 2020). Therefore, understanding the burden of the COVID-19 pandemic on mortality requires higher spatial resolution and models that account for spatial, temporal and spatio-temporal dependencies.
When working at high spatio-temporal resolution, data are generally sparse, leading to excess mortality estimates that are highly variable. This is aggravated by the fact that they are subject to spatial and temporal correlation, making it essential to use statistical methods that account for these dependencies in order to provide robust and accurate estimates. The disease mapping framework, which is commonly employed in epidemiology to study the spatio-temporal variation of diseases (Waller et al., 1997;Moraga, 2018), can be exploited to estimate excess mortality at subnational and weekly level. The Bayesian approach is naturally suited in this context, as it is able to model complex dependency structures, as well as to incorporate uncertainty in the data and modelling process. In addition, while fully propagating the uncertainty, it allows us to summarise the results at any desired spatio-temporal aggregation (using for instance coarser geographical units suitable for policy implementation). This in combination with fast approximate methods to inference, such as the Integrated Laplace Approximation (INLA, (Rue et al., 2009)), make this framework particularly powerful and appealing for monitoring of the pandemic burden and timely policy making.
In this tutorial, we describe how to run a study on excess mortality at high spatial and temporal resolution using a Bayesian approach and R. This tutorial provides a detailed implementation of the approach followed previously in 5 European regions (Konstantinoudis et al., 2022) with a slightly modified main model (Riou et al., 2023). Figure 1 illustrates graphically the workflow followed in this paper together with the main r-packages used. Source code for replicating the data wrangling (R-files starting with 01, 02 or 03), analysis (R-files starting with 04) and post-processing steps (Rfiles starting with 05 or 06 and the App folder) and data files are available from GitHub at https: //github.com/gkonstantinoudis/TutorialExcess. This tutorial is structured as follows: we first describe the modelling framework and present the case study in Italy. We then show how to run and evaluate the model, and extract and visualise the results. Finally, we present an R-shiny app which makes the results effectively and easily disseminated.

Bayesian hierarchical spatio-temporal model to estimate excess mortality
We propose a Bayesian hierarchical model to quantify spatio-temporal excess mortality under extreme events such as the COVID-19 pandemic, stratified by specific age-sex population groups. To do so, we first describe the statistical model for predicting the number of deaths from all-causes based on historical data, in the counterfactual scenario in which the pandemic did not take place. Then, we show how to estimate the magnitude of excess deaths over a specific period of time, with associated uncertainty, by comparing the predicted versus the actual number of deaths.
Let y jstk and P jstk be the number of all-cause deaths and the population at risk for the j-th week, j = 1, . . . , J t , where J t is the total number of weeks of the year t (t = 2015, . . . 2019), the s-th spatial unit (s = 1, . . . S, where S is the number of provinces in Italy), and k-th age-sex group (k = 1, . . . 10). Also let x 1jt , x 2t and z jst denote the adjustment covariates, respectively public holidays (i.e., x j = 1 if week j contains a public holiday and 0 otherwise), the year that the j-the week belongs to and temperature. We assume that the historical observed number of deaths, conditional on the risk r jstk , follows a Poisson distribution, with a log-linear model on the risk. To simplify notation, we omit k, although the following model was fitted to all k age-sex groups: y jst |r jst ∼Poisson(µ jst = r jst P jst ), Here, β 0jt is the week specific intercept in year t given by β 0jt = β 0 + jt for the k-th age-sex group, where β 0 is the global intercept and jt is an unstructured random effect representing the deviation of each week from the global intercept, which is modelled using independent and identically (iid) distributed Gaussian prior distribution with zero-mean and variance equal to τ −1 . The parameter β 1 and β 2 are unknown regression coefficients associated to the public holidays covariate x 1jt and a long term linear trend. The effect of temperature is allowed to be non-linear f (·) by specifying a second-order random walk (RW2) model: Terms b s and w j are spatial and temporal random effects, respectively. We specify the spatial random effect term, b s , with a convolution prior (Besag et al., 1991), and the temporal random effect term, w j with a non-stationary in time prior. In detail, we model b using a reparameterisation of the popular Besag-York-Mollié prior, which is a convolution of an intrinsic conditional autoregressive (CAR) model and an iid Gaussian model. Let u s be the spatially structured component defined by an intrinsic CAR (Besag, 1974) whereū is the local mean and ∂ s and m s are respectively the set and the number of neighbours of area s, τ −1 u the conditional variance, and v s the unstructured component with prior v s iid ∼ Normal(0, τ −1 v ). We model b s as follows (Besag et al., 1991;Riebler et al., 2016;Konstantinoudis et al., 2020): where u s and v s are standardised versions of u s and v s to have variance equal to 1 (Simpson et al., 2017). The term 0 ≤ φ ≤ 1 is a mixing parameter, which measures the proportion of the marginal variance explained by the structured effect. Finally, we assign to the temporal random effect term, w jt , a Gaussian random walk model of order 1 (RW1). This component captures seasonality and is specified as: The Bayesian representation of the above model is completed once we select priors for the fixed effects β 0 and β and the hyperparameters: τ , τ z , τ b , τ w , and φ. For the fixed effects we selected minimally informative Normal distributions whereas, we specified "penalising complexity" (PC) priors (Simpson et al., 2017) for the hyperparameters. PC priors are defined by penalising deviations from a "base" model (e.g., specified in terms of a specific value of the hyperparameters) and have the effect of regularising inference, while not implying too strong information. Technically, PC priors imply an exponential distribution on a function of the Kullback-Leibler divergence between the base model and an alternative model in which the relevant parameter is unrestricted. This translates to a suitable "minimally informative", regularising prior on the natural scale of the parameter.
In order to quantify the weekly excess mortality at sub-national level for specific age-sex population groups, we need to predict the number of deaths had COVID-19 not occurred. In Bayesian analysis, this can be performed by drawing random samples from the posterior predictive distribution (that is, the distribution of unobserved values conditional on the observed values from previous years). Specifically, letting θ θ θ be the model parameters, D be the observed data, and y jst * be the count of deaths that we want to predict, we have: Operationally, we first generate random samples from the joint posterior marginal of the fitted linear predictor specified in equations (1) at the highest spatial resolution available (NUTS3 regions; Nomenclature of Territorial Units for Statistics 3 regions, https://ec.europa.eu/eurostat/web/nuts/ background/). Successively, we use these to be in turn the mean parameter of a Poisson distribution, to obtain the predicted number of deaths, which represents the baseline number of deaths assuming the pandemic did not take place.
Finally, to estimate the magnitude of excess deaths, the predicted number of deaths is compared against the actual observed number of deaths, allowing the computation of the relative change in the mortality (i.e., relative to what we could expect if the pandemic did not occur). This is obtained by (i) subtracting the predicted number of deaths from the observed number of deaths in each time point j in the t * -th year and spatial unit s (number of excess deaths or NED), and (ii) dividing NED by the predicted number of deaths for each sample and multiplying by 100 (% relative excess mortality or REM).
Bayesian inference for the model is computed using Integrated Nested Laplace Approximation (INLA; (Rue et al., 2009), which performs approximate Bayesian inference on the class of latent Gaus-sian models (Rue and Held, 2005). Unlike simulation based Markov chain Monte Carlo method, INLA is a deterministic algorithm, which employs analytical approximations and efficient numerical integration schemes to provide accurate approximations of the posterior distributions in short computing times. The INLA software is provided through the R package INLA, which can be downloaded from https://www.r-inla.org/.

Outcome data
We retrieved all-cause mortality data during 2015-2020 in Italy from the Italian National Institute of Statistics (https://www.istat.it/). Data were available weekly (ISO week), by age (5-year age groups), sex and NUTS3 regions. As the COVID-19 mortality rates increase with age, we aggregated mortality counts based on the following age groups: <40, 40-59, 60-69, 70-79 and 80 years and older (Davies et al., 2021).

Population data
Population data in Italy during 2015-2020 were retrieved from the Italian National Institute of Statistics. The data represent the population in Italy on January 1st of every year and are available by age (5-year age groups), sex and NUTS3 regions. To retrieve weekly population, we performed linear interpolation by the selected age groups (<40, 40-59, 60-69, 70-79 and 80+), sex and NUTS3 regions using populations at January 1st of the current and the next year. Population counts on January 1st 2021, which takes COVID-19 deaths in 2020 into consideration, were available at the time of analysis. Our goal was, however, to predict mortality for 2020 had the pandemic not occurred. Thus we performed an additional linear interpolation by age, sex and NUTS3 regions to predict the population at January 1st 2021, using the years 2015-2020 (Figure 2, panel A). Object pop is a tibble containing the NUTS3 region ID (NUTS318CD), age group (ageg), sex (sex), year (year) and population counts ( We can use the following code (based on tidyverse and "piping" principles) to calculate the population of the 1st of January 2021 by NUTS3 regions, sex and age: pop %>% group_by(NUTS3, sex, age) %>% summarise(pop = as.vector(coef(lm(pop~year)) %*% c(1, 2021))) %>% mutate(year = 2021) -> pop2021 We acknowledge that the linear trend in the population is a rather simplistic assumption. In subsequent analyses in Switzerland, we proposed a spatio-temporal approach similar with (1) to model the population counts had the pandemic not occurred (Riou et al., 2023). The code for that analysis is also online available online (https://github.com/jriou/covid19_ascertain_deaths). Once we obtained the year 2021 we performed an additional linear interpolation to calculate weekly number of population as shown on Figure 2, panel B.

Covariates data
We used covariates related with ambient temperature and national holidays and year of death to help the model predictions. Data on air-temperature during 2015-2020 in Italy at 2m above the surface of land were retrieved from the ERA5 reanalysis data set of the Copernicus climate change program (Hersbach et al., 2020). The geographical resolution of the ERA5 estimates is 0.25 • × 0.25 • (panel A of Figure 3). We calculated the weekly mean by the centroids of the 0.25 • × 0.25 • grid (panel B of Figure  3) and then averaged the weekly temperature over the ERA5 centroids that overlay with the NUTS3 regions (panels B and C of Figure 3).

Fitting the model
The modelling process in INLA consists of three main steps: (1) the selection of priors; (2) definition of the model "formula" (which sets out the expression for the generalised linear predictor); and (3) the call to the main function inla, which computes the estimates.
In particular, we constructed the PC priors for σ = √ 1/τ , σ z = √ 1/τ z , σ b = 1/τ b and σ w = √ 1/τ w based on statement that it is unlikely to have a relative risks higher than exp(2) based solely on spatial, yearly and seasonal variation, Figure 4, panel A. For the mixing parameter φ, we set  Pr(φ < 0.5) = 0.5 reflecting our lack of knowledge about whether overdispersion or strong spatial autocorrelation should dominate the field b, Figure 4, panel B.

Model validation
To examine model validity, we performed a cross validation leaving out one historical year at a time and predicting the weekly number of deaths by NUTS3 regions for the year left out. For each stratum, we calculated the correlation between observed and fitted and a coverage probability, i.e. the probability that the observed death fall into the 95% credible interval (95% CrI) of the predicted.

Cross validation
We overall found good predictive ability of our model.
Panels 1B and 1C of Figure 5 show the median temporal nationwide excess together with 95% CrI by sex after aggregating the different age groups (using d_week and the "sex" selection). We observe a clear first pandemic wave during March and May and a second one during mid October and December in 2020. During the first pandemic wave, there were weeks when the median REM reached almost 100% in males, Figure 5. Panels 2B, 2C, 3B and 3C of Figure 5 show the median spatio-temporal excess together with 95% CrI by sex after aggregating over the different age groups. Panels 2B and 2C highlight the region of Puglia, where during the first wave of the pandemic experienced a positive, similar with the nationwide excess. When we increase the spatial resolution in panels 3B and 3C, we highlight the province of Foggia, where there was insufficient evidence of a positive excess during the first wave, but strong during the second.

Shiny Web-Application
To be able to effectively examine and communicate the different aggregation levels of the output of our modelling framework, we have also developed a Shiny Web-Application (WebApp), Figure 7. The WebApp provides spatial, temporal and spatio-temporal analysis tabs, and within each tab there are plots and summary statistics for the level of aggregation selected from the drop-down menu. Users can select across different variables (REM or NED), statistics (median or posterior probability), sex (males, females or both), age group (40 <, 40 − 59, 60 − 69, 70 − 79, 80 > and all) and different geographical level (national, NUTS2 or NUTS3 regions). Summary statistics for each area are available and they are displayed in a pop-up window, which is activated by clicking on the area of interest. In addition, graphical pop-ups are provided to show the weekly estimates for each area with the leafpop R-package (Appelhans and Detsch, 2021), in the spatio-temporal analysis tab. The WebApp that we have developed is hosted at http://atlasmortalidad.uclm.es/italyexcess/.

Summary
This tutorial provides a detailed implementation of the framework followed previously to calculate excess mortality during the COVID-19 pandemic in 5 European regions (Konstantinoudis et al., 2022). The main model used here is slightly modified based on updated results (Riou et al., 2023). We have proposed a Bayesian framework for estimating excess mortality and shown how to use R and INLA to retrieve fast and accurate estimates of the excess mortality. The proposed framework also allows combining different models and presenting the results in any age, sex, spatial and temporal aggregation desired. We have given a practical example of how to use the proposed framework modelling the excess mortality during the 2020 COVID-19 pandemic in Italy at small area level. We also developed a Shiny App to effectively communicate the results. This framework can be extended to monitor mortality from other extreme events for instance natural hazards such as hurricanes (Acosta and Irizarry, 2020). Other extensions include different ways of modelling the younger age groups to increase the predictive ability, for instance using a zero-inflated Poisson. All the above make the proposed framework particularly flexible, powerful, generalisable and appealing for online monitoring of the pandemic burden and timely policy making.